Solving the m-mixing problem for the three-dimensional time-dependent Schrodinger 
equation by rotations: application to strong-field ionization of H 2 + 
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We present a very efficient technique for solving the three-dimensional time-dependent Schrodinger 
equation. Our method is applicable to a wide range of problems where a fullly three-dimensional 
solution is required, i.e., to cases where no symmetries exist that reduce the dimensionally of the 
problem. Examples include arbitrarily oriented molecules in external fields and atoms interacting 
with elliptically polarized light. We demonstrate that even in such cases, the three-dimensional 
problem can be decomposed exactly into two two-dimensional problems at the cost of introducing 
a trivial rotation transformation. We supplement the theoretical framework with numerical results 
on strong-field ionization of arbitrarily oriented H2 4 " molecules. 
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I. INTRODUCTION 

In atomic physics the spherical symmetry of atoms 
promotes the spherical coordinates to a special position. 
The three independent variables are (r,9,<fi), with r the 
radial distance of the electron with respect to the nu- 
cleus, 9 the polar angle and 4> the azimuthal angle. The 
Schrodinger equation for the hydrogen atom is separa- 
ble in these coordinates with wave functions of the form 
ip n i m { r ) = Rni(r)Yi m (6,<fi) separated into a radial wave 
function R n i(r) and a spherical harmonic Y/ m (#, cf>). Fur- 
thermore, configurations of this type with associated or- 
bitals ipnim{r) form the building blocks of Slater deter- 
minants and consequently of mean field approaches to 
atomic structure. Even for molecules where the pres- 
ence of multiple nuclei breaks the spherical symmetry 
([L 2 ,.ff] 7^ 0), single-centre expansions in spherical har- 
monic basis has been used successfully [l|. 

For a general problem involving a single active elec- 
tron we are thus led to the consideration of the three- 
dimensional time-dependent Schrodinger equation in 
spherical coordinates and for the reduced wave function 
($ = r^>) we seek a solution of the form 



(1) 



1=0 : 



A very important advantage of this representation is 
that we can benefit from angular momentum theory 
when dealing with the angular degrees of freedom. An 
outstanding problem, however, remains. The problem, 
which is referred to as the m-mixing problem among com- 
putational scientists, is that often couplings — external 
or internal — are present that introduce a mixing of m's 
across Vs. Such m-mixings occur for example when an 
atom is subject to an elliptically polarized field or to a lin- 
early polarized field described beyond the dipole approx- 
imation. When m is no longer conserved, the dynamics 
affects all three coordinates and a numerical simulation 
is difficult: three-dimensional calculations tend to be ex- 



tremely time-consuming and computationally demand- 
ing. 

In the course of our recent work concerned with 
alignment-dependent response of molecules to strong ex- 
ternal fields we found a solution that speeds up the cal- 
culation by the use of an exact mapping of the three- 
dimensional problem to two two-dimensional problems. 
In the following we discuss the method by the specific 
example of the response of an arbitrarily oriented di- 
atomic molecule to an external perturbation so strong 
that the system is ionized. As will become clear, the cen- 
tral ideas are completely general and carry over to the 
related case of atoms in elliptically polarized fields, poly- 
atomic molecules as well as m-problems in geology and 
astronomy where expansions in spherical harmonics are 
also often encountered. 

The paper is organized as follows: in Sec. |TT] we give 
an overview over the basic idea of our technique. In 
Sec. Mil we outline the numerical implementation and 
discuss physical results for H2 + strong-field ionization. 
Sec. IIVI concludes. 



II. BASIC IDEAS AND PRINCIPLES 

We illustrate the basic ideas and principles of the 
method by discussing the specific example of a linear di- 
atomic molecules in an external electromagnetic field. In 
Fig. [1] we show the coordinate systems which are rel- 
evant for the field-molecule problem. The coordinates 
(x L ,y L , z L ) specify the laboratory (L) fixed coordinate 
system defined by the external polarization vector. We 
assume that the field is linearly polarized and return to 
the generalization to elliptically polarized light in Sec. IIVI 
The coordinate system denoted by superscripts M is the 
molecular fixed frame and is rotated by the Euler angles 
(a, (3, 7) with respect to the laboratory fixed system. The 
rotation is accomplished by an a rotation around the z L - 
axis, followed by a (3 rotation around the y M -axis, and 
finally a 7 rotation around the z M -axis. For the case 
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FIG. 1: The orientation of the molecular coordinate system 
(M, dashed) with respect to the laboratory fixed system (L, 
solid). In the figure only the Euler angles a, /3 are nonvanish- 
ing. 



considered the only really distinct geometries are associ- 
ated with the angle f3. Results for different orientations 
due to the angle a are trivially related by a simple rota- 
tion around the z L axis. Also the 7 rotations around the 
molecular axis are insignificant as a consequence of the 
axial symmetry of the molecule. 

We want to determine how the wave function of an 
electron is affected by the operators and V^(t), 

corresponding to the interaction with the nuclei and the 
field, respectively. We assume that we can treat these two 
operators separately, which is the case in a split-operator 
approach as described in Sec. [TTT] below. Our strategy 
is first to represent the wave function in the molecular 
frame and calculate the action of V^ M \ Secondly, we 
transform the updated wave function to the laboratory 
fixed frame and apply the operator V^ 1 '. Finally we can 
return to the molecular frame by the inverse rotation. 
These forward (/3) and backward (— (3) rotations of the 
wave function are illustrated in Fig. [21 The active inter- 
action (y( M > or V^) is marked by black and the inactive 
operation is gray. This propagation scheme for arbitrary 
orientation of the polarization axis with respect to the 
internuclear axis, exhibits the strength of the present ap- 
proach since it allows us to perform the calculations very 
efficiently. Whenever we apply an axially symmetric op- 
erator, we do not mix different m states provided that the 
wave function is expressed in the proper reference frame. 
Thus we can apply the operator separately on each differ- 
ent m state. The decoupling of different m states means 
effectively that we have reduced the three-dimensional 
problem to a number of two-dimensional problems in ad- 
dition to two rotation operations. 

The rotation transformation is in principle possible in 
all sets of coordinates and the separation in m applies 
to any coordinate system where the azimuthal angle 4> 
is an independent variable, e.g. cylindrical, parabolic, or 
spheroidal coordinates. The two unique features of the 



FIG. 2: Schematic picture of the rotation operation. The 
contour lines indicate the field free 1<t 9 ground state of H2 4 " 
in the xz plane. The double headed arrow shows the direction 
of the laser polarization vector. In (a) we calculate the action 
of the molecular potential and express the wave function in 
the molecular frame with the internuclear axis parallel to the 
z axis. In panel (b) we transform the wave function to the 
laboratory fixed system with the laser polarization parallel to 
the z axis in order to propagate by the field interaction. The 
transformation between the two frames is represented by the 
rotation operator D. 



spherical representation JT]) are that (i) the transforma- 
tion matrix contains Wigner rotation functions which are 
known analytically and (ii) the transformation is guaran- 
teed to be exactly unitary for functions that are band- 
width limited by a maximum I — Z max , i.e., the popula- 
tion in states with I > Z max is zero. 



III. NUMERICAL RESULTS 

In the present work, we solve the time-dependent 
Schrodinger equation (TDSE) for the electronic motion 
in H2 + in the presence of a time-dependent electromag- 
netic field. We represent the angular variables in a basis 
of spherical harmonics and write the reduced wave func- 
tion as in Eq. JT]). The radial functions fi m which contain 
the time dependence are discretized on an equidistant 
spatial mesh. The expansion in spherical harmonics is 
truncated such that I < / max leading to a total number 
of (/max + l) 2 angular basis functions. The reduced wave 
function satisfies the TDSE with the Hamiltonian [atomic 
units (h = \e\ = m e = ao = 1) are used throughout] 

H ( t ) = -\^ + ^ + V(t)=T r + T l + V(t), (2) 

where L is the usual angular momentum operator and V 
includes the electronic interaction with the field and the 
nuclei. We solve the time-evolution from time t to t + r 
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numerically by using the split-operator technique 
$(r,i+r) = e-^e-^e-^^+^e-^ie-^^r.t). 

The error in the propagation scheme above is approxi- 
mately cubic in r and occurs mainly due to the split- 
ting of non-commuting operators. A related propagation 
scheme was applied in geometries with azimuthal sym- 
metry 0], and the propagation techniques used for the 
kinetic operators T r and Tj are readily extended to our 
fully three-dimensional problem. We will therefore turn 
to the new propagation method of the molecular poten- 
tial and the field interaction. 

We describe the electromagnetic field in the dipole ap- 
proximation by the vector potential 

A(t) = eA a (t)cos{ujt), (4) 

where Ao(t) is the envelope function, ui the frequency and 



e the polarization direction. The electric field is obtained 
as F(t) = -dA(t)/dt, The operator V in Eq. © is writ- 
ten as the sum of the field interaction and the molecular 
potential 

V(t) = V$ L (t)+vX, (5) 

where the subscripts denote the variables on which the 
operators act. Om is the polar angle in the molecular 
frame [Fig. [2] (a)] and 9l the polar angle in the labora- 
tory fixed system [Fig. [2] (b)]. The molecular operator is 
diagonal in coordinate space 

^(M) =y( M) Mi)j (6) 

while the field interaction can be represented either in 
the length- (LG) or the velocity gauge (VG) as 



F(t)r cos 9 L 
iA(t) [i ( 



cos ujj + sin Ul 



do L 



COS uljt- 



LG 
VG 



(7) 



To calculate the action of V in the propagation we make 
the split 



-*v(tH 



-iV< M > T/2 e -iV< r > (t+5)r e -iV< M 'r/2 



(8) 



For each radial grid point Ti we write the wave function 
as a vector in the spherical harmonics basis, cf. Eq. ^ 



f {M) {n,t) 



ATM 



V 



(9) 



where the coefficients refer to the molecular frame. The 
molecular potential is diagonal in the radial coordinate, 
and cannot induce mixings vectors that belong to dif- 
ferent radial coordinates. We evaluate the action of 



-iV< M V/2 



by its matrix representation in the spherical 



harmonics basis for each fixed value of r 

(lm\e- iV ' M) ^ 2 \l'm') = 5 mm ,{lm\e- ivW ^ e ^/*\l'm). 

< (10) 

The selection rule m = m! occurs since ' is indepen- 
dent of 4>m- Now it is evident that e~ iy( >r / 2 is repre- 



sented by a the block diagonal form 



m=0 
1=0,1,2. ■■■ ,l„ 







1=1,2, 



\ 







[ 



l=l„ 



V 



(11) 

Although not essential for our present discussion, we note 
that for inversion symmetric potentials as in the case of 
H2 + , a further block diagonalization in even and odd 
parity blocks can be obtained. From the block diagonal 
structure of the matrix representation, it is clear that the 
propagation can be accomplished separately within each 
to subspace, and the full three-dimensional propagation 
effectively reduces to independent two-dimensional prop- 
agations, which can be solved by matrix multiplications 
on each m block. There is a total number of 2Z max + 1 
individual m blocks with dimensionality between 1 and 

^max 1 • 

After having applied the molecular potential we trans- 
form the wave function to the laboratory fixed frame. 
We relate the expansion in spherical harmonics in dif- 
ferent frames by representation of the rotation operator 
in spherical harmonics, i.e., the Wigner rotation matrix 
D(a, /3, 7). The laboratory fixed expansion coefficients 
are then obtained as f {L) {n,t) = D(a,j3, 7) • / (M) (n, t). 
We note that this matrix multiplication is very fast since 
the rotation does not mix different Vs and D(a,0, 7) is 
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consequently sparse. Also note that the rotation opera- 
tion is independent of the radial coordinate and we can 
therefore use the same rotation operation on all the vec- 
tors ([H]) for different r's. 

Having obtained the wave function in the laboratory 
fixed frame, we can easily apply the field interaction op- 
erator. Again, without m-couplings, the individual two- 
dimensional problems can be solved straightforwardly [|| . 
Finally we return to the molecular frame by the inverse 
transformation / (M) (r l7 i) = 7) • / (i) (r,,i). 

We close this section with a few remarks on the scaling 
of the computations with the size of the problem. In 
an alternative three-dimensional approach where we in a 
single step treat the total V and mix between all (7 m ax + 
l) 2 angular basis states, the computational complexity 
scales as O(^max) 13 • Cur present method, on the other 
hand, scales more favourably as 0{l^ ayi ). In numerical 
simulations for typical bandwidths of Z max ~ 15 — 39, we 
have checked that both three-dimensional methods agree 
in their predictions but with a great speed-up of the order 
of a factor of 100 — 500 in favor of the new method. 



(a) (b) (c) 
1.1x10* 




FIG. 3: (Color online) Angular differential ionization proba- 
bility dP/dQ.L for the alignment angles (a) 0°, (b) 45°, and (c) 
90°. The laser polarization direction is vertical in all panels 
and the molecular axis is indicated by the thick solid line. The 
numbers on the axes indicate dP/dD.L 1^=0- The parameters 
of the electromagnetic field are specified in the text. 



A. Ionization of H2 

We calculate the ionization probability for H2 + in- 
duced by a strong infra-red light source. The two pro- 
tons are fixed at the equilibrium internuclear distance of 
2 a.u.. The field is taken to be linearly polarized with fre- 
quency u> = 0.057 a.u. (A = 800 nm), and peak intensity 
5 x 10 14 W/cm . We use a sine-square envelope func- 
tion that encloses seven optical cycles, corresponding to 
a total pulse duration of 19 fs. Convergent results are 
obtained with l max = 23 and 1024 radial grid points ex- 
tending to a box size of 150 a.u.. In order to avoid re- 
flections at the edge of the box, we impose an absorbing 
boundary. The time step size is r = 5 x 10~ 3 a.u.. We 
choose the velocity gauge form of the interaction since it 
is superior to the length gauge introducing converged 
results for dynamical problems 0, Q . 

First we calculate the angular differential ionization 
probability. For that purpose we need the gauge invariant 
current density 

J(r, t) = Re [**(r, t) (p + A(t)) #(r, t)} , (12) 

where p — — iV is the canonical momentum. We relate 
the outgoing radial probability flux at some large dis- 
tance 1Z to the differential ionization probability in the 
laboratory fixed frame 

dP f°° 

— =/ dtf L - J(K,tl L ,t)1l 2 . (13) 
uMl Jo 

We must of course choose TZ to be smaller than the ra- 
dial distance at which we turn on the absorbing potential. 
Figure [3] shows the angular differential probabilities for 
the alignment angles 0°, 45°, and 90°. In all cases, the 
electron escapes exclusively in a very narrow cone along 



the polarization direction. These results are in accor- 
dance with expectations from the quasistatic tunneling 
picture. The ionization dynamics is often considered as 
being tunneling- like for strong, low frequency fields where 
the Keldysh parameter fullfils 7 < 1 [7]. In the present 
case 7 = 0.7 at the peak intensity. In the tunneling pic- 
ture the electron is assumed to escape near the field di- 
rection since the barrier has its shortest spatial extension 
in that direction Q . 

The most notable difference between panels (a)-(c) is 
the overall scaling of the distribution which decreases 
with increasing angle between the polarization and in- 
ternuclear axes. We can qualitatively explain this ob- 
servation by the associated decrease in electronic charge 
density of the intial <7 g -orbital after the polarization di- 
rection (see countour plot in Fig. [5]). The same reason- 
ing carries over to the behavior of the total alignment 
dependent ionization probabilities shown in Fig. [U The 
results in this figure can be obtained by integrating the 
differential ionization probability Eq. (|13[) over all direc- 
tions. Alternatively, we may project out the bound state 
components of the final wave function. For comparison, 
Fig. Q] also contains the results from Ref. @ which were 
obtained by a field of the same frequency and peak inten- 
sity but with a slightly different pulse shape (trapezoidal) 
and longer duration. We find somewhat lower ionization 
probabilities than in Ref. |6| since our pulse is at the peak 
intensity for a shorter duration of time. Although the 
two data series are not directly comparable, the overall 
behaviour is similar, namely decreasing ionization prob- 
ability with increasing alignment angle from parallel (0°) 
to perpendicular (90°). 
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FIG. 4: Total ionization probability as a function of alignment 
angle. The present results are indicated by the solid line. The 
dashed line is taken from Ref. |6j after scaling by the factor 
0.18. 



IV. CONCLUSION AND OUTLOOK 

In conclusion, we have developed a new approach that 
accurately and efficiently resolves the m-mixing prob- 
lem in large scale computations in a spherical coordinate 
system. The method relies on an identification of ro- 
tations in the intermediate propagation that brings the 
wave function into a frame of reference in which m is 
conserved. This means that time-consuming m-mixing 
induced by the external perturbation is avoided and in- 
stead delegated to the rotations which are very efficiently 



implemented using the Wigner rotation matrix represen- 
tation of the rotation operator in the spherical harmonics 
basis. 

We have chosen the linear molecule interacting with 
a linearly polarized field to illustrate our method, but a 
similar approach can be used in a much broader range of 
three-dimensional problems. For example we could con- 
sider an elliptically polarized field. In the split operator 
method we take the time step r to be small enough such 
that the field can be taken to be constant both in mag- 
nitude and polarization direction within the small time 
interval. We can therefore consider a time-dependent 
laboratory frame which follows the instantaneous polar- 
ization direction. If we make the transformation from 
the molecular frame to the new laboratory frame, we are 
again able to treat the field as being linear and propagate 
as discussed above. Our method can also be extended to 
arbitrary nuclear positions. For any nuclear configura- 
tion, we can attach a coordinate system to each nucleus 
with a z axis from the origin to the nucleus. Then we de- 
compose the molecular potential to a sum of nuclear po- 
tentials, each of which can be propagated with azimuthal 
symmetry in their own reference frame. Despite the fact 
that we now need rotations between the coordinate sys- 
tems belonging to all of the nuclei, the total calculation is 
still in the same complexity class with respect to scaling 

in /max- 
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